function [N] = ElecDens(s, perigee, data, lat, lon)
%
% Electron Density
%
%DESCRIPTION:
%This function implements the computation of the Electron Density, sorting,
%on the base of the altitude, if the bottom or the top side N should be
%considered.
%
%PROTOTYPE:
% [N] = ElecDens(s, perigee, data)
% [N] = ElecDens(s, perigee, data, lat, lon)   (see NOTES)
%
%--------------------------------------------------------------------------
% INPUTS:
%   s          [1x1]       Point of Computation      [km]
%   perigee    [---]       Ray-Perigee Data          [strc] (see NOTES)
%   data       [---]       Problem Data              [strc] (see NOTES)
%   lat        [1x1]       Latitude                  [deg] (optional)
%   lon        [1x1]       Longitude                 [deg] (optional)
%--------------------------------------------------------------------------
% OUTPUTS:
%   TEC        [1x1]       Total Electron Content    [TECU]
%--------------------------------------------------------------------------
%
%NOTES:
% - For the tollerance (accuracy) "tol" the adviced value is 0.01.
% - The input "perigee" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       perigee.rp      = Radius           [km]
%       perigee.latp    = Latitude         [deg]
%       perigee.lonp    = Longitude        [deg]
%       perigee.sinlatp = Sine of Latitude [-]
%       perigee.coslatp = Cosine of Latit. [-]
%       perigee.sinsigp = Sine of Zenith   [-]
%       perigee.cossigp = Cosine of Zenith [-]
% - The input "data" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       data.a0      = 1st Az Coeff   [-]
%       data.a1      = 2nd Az Coeff   [-]
%       data.a2      = 3rd Az Coeff   [-]
%       data.mth     = Month          [month]
%       data.UT      = Universal Time [hours]
%       data.stModip = MODIP Table    [-] (modipneqg_wrapped.asc)
%       data.F2      = F2 Table       [-] (ccir21.asc)
%       data.Fm3     = Fm3 Table      [-] (ccir21.asc)
% - The inputs "lat" and "lon" are optional and shall be given ONLY for
%   Vertical TEC computation!!!
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Great Circles
if nargin == 3
    [hs, lats, lons]  = CoordPath(s, perigee);
elseif nargin == 5
    hs = sqrt(s^2 + perigee.rp^2) - 6371.2;
    lats = lat;
    lons = lon;
else
    disp('Uncorrect number of inputs! Check ElecDens.m description.');
    return
end

%Ionospheric Model
[mu]                = MODIPcomp  (lats, lons, data.stModip);
[Az]                = EffAz      (data.a0, data.a1, data.a2, mu);
[Az_R]              = EffSSN     (Az);
[AF2, Am3]          = interpAz_R (data.F2, data.Fm3, Az_R);
[CF2, Cm3]          = FourierF2  (data.UT, AF2, Am3);
[foF2, M3000, NmF2] = LegendreF2 (mu, lats, lons, CF2, Cm3);
[sinSun, cosSun]    = SolDecl    (data.mth, data.UT);
[LT]                = LocTime    (data.UT, lons);
[chi]               = SolZenith  (lats, LT, sinSun, cosSun);
[chi_eff]           = EffSolZen  (chi);
[foE, NmE]          = foE_NmE    (lats, Az, chi_eff, data.mth);
[hmF2]              = hMF2       (foE, foF2, M3000);
if     hs <= hmF2
    %Bottomside N
    [hmE]        = hME;
    [hmF1]       = hMF1(hmF2, hmE);
    [A1]         = A_1(NmF2);
    [foF1, NmF1] = foF1_NmF1(foE, foF2);
    [B2bot]      = B2BOT(NmF2, foF2, M3000);
    [B1top]      = B1TOP(hmF1, hmF2);
    [B1bot]      = B1BOT(hmF1, hmE);
    [BEtop]      = BETOP(B1bot);
    [BEbot]      = BEBOT;
    [A2, A3]     = A_23(NmE, NmF1, A1, hmF2, hmF1, hmE, BEtop,...
                       B1bot, B2bot, foF1);
    
    [N] = TEC_BOT(hs, A1, A2, A3, hmF2, hmF1, hmE, B2bot, B1top,...
                                      B1bot, BEtop, BEbot);
elseif hs > hmF2
    %Topside N
    [B2bot] = B2BOT(NmF2, foF2, M3000);
    [k]     = ShapeK(data.mth, NmF2, hmF2, B2bot, Az_R);
    [H0]    = Height0(B2bot, k);
    
    [N] = TEC_TOP(hs, NmF2, hmF2, H0);
end

end